Predicción de rendimientos del S&P500¶

- Modelo SkForecaster con Random Forest-¶

In [1]:
#!pip install skforecast

#Librerías para cálculos y gráficos
#====================================================================================
import pandas as pd
import numpy as np
from numpy import array
import plotly.graph_objs as go
import matplotlib.pyplot as plt
#pip install matplotlib==3.7.1
import seaborn as sns

#Librerías para Preprocesamiento de datos
#====================================================================================
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer


#Librerías para Implementación de Modelos
#====================================================================================
from sklearn.ensemble import RandomForestRegressor
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregCustom import ForecasterAutoregCustom
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import backtesting_forecaster
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from skforecast.utils import save_forecaster
from skforecast.utils import load_forecaster

Importamos la base de datos con la cual vamos a trabajar

In [2]:
dfmacro = pd.read_csv(r'C:\Users\Vivi\Downloads\BaseAccionesF.csv')
dfmacro.columns = ["Date","GSPC","AAPL","AMZN","MSTF","TSLA","GOOG","GOOGL","NVDA","BRK.B","META","UNH","JNJ","PG","VIX","UNRATE","FEDFUNDS","UMCSENT","CPIAUCSL","DOLARINDEX"]
dftotal = pd.read_csv(r'C:\Users\Vivi\Downloads\BaseAcciones2.csv')
dftotal['Date'] = dfmacro['Date']
dftotal = dftotal.dropna()
In [3]:
dftotal.iloc[:,[43,56,75]].describe()
Out[3]:
GSPC.ReturnSuavizado VIX.ReturnSuavizado DolarIndex.ReturnSuavizado
count 2667.000000 2667.000000 2667.000000
mean 0.062439 -0.344104 0.008501
std 0.880420 7.081882 0.298783
min -3.002265 -26.622753 -1.188374
25% -0.360074 -4.487971 -0.166130
50% 0.063567 -0.749383 -0.001720
75% 0.544763 3.449533 0.181680
max 3.101521 26.760500 1.222749
In [4]:
df2= dftotal.loc[:,['Date','GSPC.ReturnSuavizado', 'AAPL.ReturnSuavizado', 'AMZN.ReturnSuavizado',
       'MSTF.ReturnSuavizado', 'TSLA.ReturnSuavizado', 'GOOG.ReturnSuavizado',
       'GOOGL.ReturnSuavizado', 'NVDA.ReturnSuavizado',
       'BRK.B.ReturnSuavizado', 'META.ReturnSuavizado', 'UNH.ReturnSuavizado',
       'JNJ.ReturnSuavizado', 'PG.ReturnSuavizado', 'VIX.ReturnSuavizado','DolarIndex.ReturnSuavizado']]
df2['Date'] = pd.to_datetime(df2.Date)
df2 = df2.set_index('Date')
df2 = df2.asfreq('B',method='bfill') 
#df2 = df2.iloc[:,1:]
print(f'Número de filas con missing values: {df2.isnull().any(axis=1).mean()}')

Número de filas con missing values: 0.0
In [5]:
# Verificar que un índice temporal está completo
# ==============================================================================
(df2.index == pd.date_range(
                    start = df2.index.min(),
                    end   = df2.index.max(),
                    freq  = df2.index.freq)
).all()
Out[5]:
True
In [6]:
# Verificar que no existan valores nulos
# ==============================================================================
print(f'Número de filas con missing values: {df2.isnull().any(axis=1).mean()}')

Número de filas con missing values: 0.0

Construcción de conjuntos de entrenamiento y prueba¶

In [7]:
# Definir tamaños de train y test
# ==============================================================================

dataindex2= df2.index

steps = 20
n100 = len(df2)
ntrainall = int(n100-steps)
ntestall = n100 - ntrainall
ntrain = int(ntrainall-steps)
ntest = steps

n100, ntrainall, ntestall,ntrain,ntest
Out[7]:
(2812, 2792, 20, 2772, 20)
In [8]:
# Construir conjuntos de entrenamiento y validación
# ==============================================================================
dftrainall = df2[:ntrainall]
dftestall = df2[ntrainall:]

len(dftrainall),len(dftestall)
Out[8]:
(2792, 20)
In [9]:
# Sacamos una porción de datos de train para validación nuevamente
# ==============================================================================
dftrain_sinest = dftrainall[:ntrain]
dftest_sinest = dftrainall[ntrain:]
dftestall_sinest = df2[ntrainall:]

len(dftrain_sinest),len(dftest_sinest),len(dftestall_sinest)
Out[9]:
(2772, 20, 20)

Estandarización de variables¶

In [10]:
scaler = StandardScaler()
In [11]:
cols_to_standardize = list(dftrain_sinest.columns)


dftrain = pd.DataFrame(scaler.fit_transform(dftrain_sinest[cols_to_standardize]),
                          index=dftrain_sinest.index, columns=cols_to_standardize)

# Estandarizamos los datos de test
dftest = pd.DataFrame(scaler.transform(dftest_sinest[cols_to_standardize]),
                              index=dftest_sinest.index, columns=cols_to_standardize)

dftestall = pd.DataFrame(scaler.transform(dftestall_sinest[cols_to_standardize]),
                              index=dftestall_sinest.index, columns=cols_to_standardize)
In [12]:
#dftrain.to_excel('dftrain.xlsx')
#dftest.to_excel('dftest.xlsx')
#dftestall.to_excel('dftestall.xlsx')

Experimento 1: Forecaster Autoregresor con variables exógenas¶

En primera instancia se corre el modelo con el 100% de las variables con las que contamos, dejando a y como predictor y en variables exógenas los rendimientos de las acciones y las variables macroeconómicas.

Se calculan los hiperparámetros y se le pide al modelo que retorne el modelo mejor entrenado para las predicciones:

In [13]:
# Cálculo de hiperparámetros por grid search
# ==============================================================================

forecaster1 = ForecasterAutoreg(
                regressor = RandomForestRegressor(random_state=123),
                lags      = 14 
             )

# Lags used as predictors
lags_grid = [5, 20]

# Regressor's hyperparameters
param_grid = {'n_estimators': [100, 500],
              'max_depth': [3, 5, 10]}

results_grid = grid_search_forecaster(
                        forecaster         = forecaster1,
                        y                  = dftrain.iloc[:,0],
                        exog               = dftrain.iloc[:,1:],
                        param_grid         = param_grid,
                        lags_grid          = lags_grid,
                        steps              = steps,
                        refit              = True,
                        metric             = 'mean_squared_error',
                        initial_train_size = int(len(dftrain)*0.5),
                        fixed_train_size   = False,
                        return_best        = True,
                        verbose            = False
                        )

Number of models compared: 12.

loop lags_grid:   0%|                                               | 0/2 [00:00<?, ?it/s]
loop param_grid:   0%|                                              | 0/6 [00:00<?, ?it/s]C:\Users\Vivi\tf\lib\site-packages\skforecast\model_selection\model_selection.py:370: RuntimeWarning: The forecaster will be fit 70 times. This can take substantial amounts of time. If not feasible, try with `refit = False`. 

  warnings.warn(

loop param_grid:  17%|██████▎                               | 1/6 [01:27<07:18, 87.69s/it]
loop param_grid:  33%|████████████▎                        | 2/6 [08:36<19:13, 288.40s/it]
loop param_grid:  50%|██████████████████▌                  | 3/6 [10:45<10:46, 215.57s/it]
loop param_grid:  67%|████████████████████████▋            | 4/6 [21:24<12:45, 382.61s/it]
loop param_grid:  83%|██████████████████████████████▊      | 5/6 [25:06<05:24, 324.87s/it]
loop param_grid: 100%|█████████████████████████████████████| 6/6 [43:20<00:00, 586.36s/it]
loop lags_grid:  50%|██████████████████▌                  | 1/2 [43:20<43:20, 2600.63s/it]
loop param_grid:   0%|                                              | 0/6 [00:00<?, ?it/s]
loop param_grid:  17%|██████▏                              | 1/6 [02:18<11:34, 138.91s/it]
loop param_grid:  33%|████████████▎                        | 2/6 [13:50<30:55, 463.84s/it]
loop param_grid:  50%|██████████████████▌                  | 3/6 [17:24<17:30, 350.02s/it]
loop param_grid:  67%|████████████████████████▋            | 4/6 [35:14<21:08, 634.33s/it]
loop param_grid:  83%|██████████████████████████████▊      | 5/6 [41:31<09:01, 541.35s/it]
loop param_grid: 100%|███████████████████████████████████| 6/6 [1:12:48<00:00, 995.51s/it]
loop lags_grid: 100%|███████████████████████████████████| 2/2 [1:56:09<00:00, 3484.60s/it]

`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  Lags: [1 2 3 4 5] 
  Parameters: {'max_depth': 10, 'n_estimators': 500}
  Backtesting metric: 0.3251563573906958
In [14]:
# Predicciones
# ==============================================================================
predicciones = forecaster1.predict_interval(steps=steps, exog=dftest.iloc[:,1:],interval = [1, 99],n_boot   = 500)

Funciones para crear gráficos y métricas¶

In [15]:
#Podemos primeramente graficar las secciones: la serie de entrenamiento, la predicción en la sección de prueba y la sección de prueba.
def grafico_comparativo(dftrain,dftest,predicciones):
    trace1 = go.Scatter(
        x = dftrain.index,
        y = dftrain.iloc[:,0],
        mode = 'lines',
        name = 'Data'
    )
    trace2 = go.Scatter(
        x = dftest.index,
        y = dftest.iloc[:,0],
        mode = 'lines',
        name = 'Prediction'
    )
    trace3 = go.Scatter(
        x = predicciones.index,
        y = predicciones['pred'],
        mode='lines',
        name = 'Reality'
    )
    layout = go.Layout(
        title = "S&P Pronóstico",
        xaxis = {'title' : "Date"},
        yaxis = {'title' : "Close"}
    )
    fig = go.Figure(data=[trace1, trace2, trace3], layout=layout)
    fig.show()
In [16]:
grafico_comparativo(dftrain,dftest,predicciones)
In [17]:
def grafico_intervalos(y,predicciones):
    fig, ax = plt.subplots(figsize=(10, 3))
    y.plot(ax=ax, label='test')
    predicciones['pred'].plot(ax=ax, label='predicciones')
    ax.fill_between(
        predicciones.index,
        predicciones['lower_bound'],
        predicciones['upper_bound'],
        color = 'red',
        alpha = 0.2
    )
    ax.legend();
In [19]:
grafico_intervalos(dftest.iloc[:,0],predicciones)
In [20]:
def grafico(y,pred_1fin):
    fig, ax = plt.subplots(figsize=(10, 3))
    y.plot(ax=ax, label='test')
    pred_1fin['pred'].plot(ax=ax, label='predicciones')
    ax.legend();
In [56]:
def desempeño(dftest,predicciones):
    error_mse = mean_squared_error(
                    y_true = dftest.iloc[:,0],
                    y_pred = predicciones['pred']
                )

    rmse = np.sqrt (error_mse)

    print(f"EL MSE del modelo en test es de: {error_mse}")
    print(f"EL RMSE del modelo en test es de: {rmse}")

    r2 = r2_score(
        y_true = dftest.iloc[:,0],
        y_pred = predicciones['pred']
        )
    
    print(f"EL R2 del modelo en test es de: {r2}")
    
    # suponga que y_true son los valores reales y y_pred son las predicciones del modelo
    error_mae = mean_absolute_error(
                    y_true = dftest.iloc[:,0],
                    y_pred = predicciones['pred']
                )

    print(f"EL MAE del modelo en test es de: {error_mae}")
In [22]:
desempeño(dftest,predicciones,)

EL MSE del modelo en test es de: 0.13204347393516816
EL RMSE del modelo en test es de: 0.36337786660055144
EL R2 del modelo en test es de: 0.9011377899512703
EL MAE del modelo en test es de: 0.36337786660055144

Con los datos obtenidos se corre el modelo sugerido con el 100% de los datos de entrenamiento y se obtiene:

Implementación del modelo del experimiento # 1 con 100% de los datos de entrenamiento:¶

In [23]:
dftrainall_est = pd.concat([dftrain,dftest])
In [24]:
# Configuramos la función de pronóstico de forecaster con la mejor combinación obteida de hiperparámetros.
# ==============================================================================
### Aquí debe ir tu codigo, 
regressor = RandomForestRegressor(max_depth= 10, n_estimators=500 , random_state=123)
### Hasta aquí modificas

forecaster = ForecasterAutoreg(
                regressor = regressor,
                lags      = 5
             )

forecaster.fit(y = dftrainall_est.iloc[:,0], exog = dftrainall_est.iloc[:,1:])
In [25]:
# Pronóstico, en este caso elegimos una ventana de 20 días.
# ==============================================================================
pred_1fin = forecaster.predict_interval(steps=steps, exog=dftestall.iloc[:,1:],interval = [1, 99],n_boot   = 500)
In [57]:
grafico_intervalos(dftestall.iloc[:,0],pred_1fin)
In [58]:
desempeño(dftestall,pred_1fin)

EL MSE del modelo en test es de: 0.12163476353394123
EL RMSE del modelo en test es de: 0.3487617575565607
EL R2 del modelo en test es de: 0.9082587299062685
EL MAE del modelo en test es de: 0.25495078018323786
In [28]:
forecaster.get_feature_importance()
Out[28]:
feature importance
0 lag_1 0.008050
1 lag_2 0.005388
2 lag_3 0.005674
3 lag_4 0.006963
4 lag_5 0.006396
5 AAPL.ReturnSuavizado 0.043067
6 AMZN.ReturnSuavizado 0.023042
7 MSTF.ReturnSuavizado 0.176339
8 TSLA.ReturnSuavizado 0.009582
9 GOOG.ReturnSuavizado 0.047308
10 GOOGL.ReturnSuavizado 0.044740
11 NVDA.ReturnSuavizado 0.026820
12 BRK.B.ReturnSuavizado 0.328805
13 META.ReturnSuavizado 0.008780
14 UNH.ReturnSuavizado 0.016280
15 JNJ.ReturnSuavizado 0.013042
16 PG.ReturnSuavizado 0.014296
17 VIX.ReturnSuavizado 0.204769
18 DolarIndex.ReturnSuavizado 0.010659

Al analizar la importancia de los predictores, encontramos que el VIX, MSTF y BRK.B son los 3 predictores de mayor importancia. Dado que las variables que ingresamos como exógenas al modelo requieren ser conocidas y que esto implicaría tener un modelo predictivo por variable, vamos a explorar la forma de reducir la cantidad de variables que tenemos que predecir.

Analicemos estas dos acciones dentro del portafolio global:

In [29]:
import seaborn as sns
fig, ax = plt.subplots(figsize=(10,6))

n = len(df2.columns)
df_cor = df2.reset_index()
df_cor = df_cor.iloc[:,1:n-2]
#corr = np.corrcoef(df_cor.T)
upp_mat = np.triu(df_cor.corr())
sns.heatmap(df_cor.corr(),  vmin = -1, vmax = +1, annot = True, cmap = 'coolwarm', mask = upp_mat)
Out[29]:
<Axes: >
In [30]:
dftotal.iloc[:,1:14].describe()
Out[30]:
GSPC.Adjusted AAPL.Adjusted AMZN.Adjusted MSFT.Adjusted TSLA.Adjusted GOOG.Adjusted GOOGL.Adjusted NVDA.Adjusted BRK.B.Adjusted META.Adjusted UNH.Adjusted JNJ.Adjusted PG.Adjusted
count 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000 2667.000000
mean 2685.892010 59.416211 70.665662 113.777254 70.480455 57.459379 57.640025 62.870372 185.321669 148.262557 212.637677 110.799351 88.421621
std 912.201173 49.474965 53.398954 91.487247 100.416735 35.656729 35.156219 74.184491 65.377359 86.401787 144.348470 35.494257 31.780452
min 1278.040039 12.046193 10.411000 21.689646 1.740000 13.924059 13.990240 2.610562 78.830002 17.730000 43.228558 45.895130 42.930496
25% 1996.284973 22.743071 19.154250 39.305872 13.896667 28.364902 28.747622 4.884724 135.935005 78.385002 96.489235 81.087364 63.413342
50% 2496.840088 37.335831 50.503502 71.427315 18.148666 48.944500 49.808498 37.803696 178.570007 143.800003 183.339828 112.401558 73.754150
75% 3257.575074 86.821614 100.522503 191.692101 66.203666 73.241749 73.216999 92.120262 218.050003 190.404998 287.365905 136.825287 115.653481
max 4796.560059 180.683853 186.570496 339.075562 409.970001 150.709000 149.838501 333.350769 359.570007 382.179993 551.479736 181.108810 159.209808

El mapa de correlación nos muestra que al incluir los precios de las acciones de empresas como Microsoft (MSFT) y Berkshire Hathaway (BRK.B) en el modelo para predecir el S&P500, podríamos obtener una mejor precisión pues dado que la correlación indica una relación lineal entre dos variables, los precios de las acciones de MSFT, BRK.B y el S&P500. Como vemos, encontrmaos una correlación positiva cercana a 1 lo que puede interpretarse como que los precios de las acciones se mueven en la misma dirección que el S&P500.

De otro lado, al revisar el portafolio de acciones, observamos que el precio promedio de estas dos acciones está dentro de las acciones con los mayores precios de las seleccionadas, lo cual, vinculándolo a la forma como se calcula el S&P500 hace que estas acciones tengan una alta improtancia en el cálculo del resultado final del S&P. Si bien, hay otras acciones con mayores precios, no tienen una correlación tan alta con el índice razón por la cual no se incluirán.

Dado que el S&P500 es un índice ponderado de capitalización bursátil en donde se encuentran las 500 empresas de EEUU con gran capitalización, las empresas con mayor capitalización bursátil deben tener una mayor influencia en el resultado del índice. Adicionalmente, son empresas que pueden ser potencialmente líderes en el mercado y bien posicionadas, lo cual puede marcar tendencias en el mercado.

Por tanto, vamos a explorar diferentes posibilidadades para correr el modelo con menos variables exógenas:

Experimento 2: Forecaster Autoregresor con variables exógenas filtradas (MSTF, BRK.B y VIX)¶

In [31]:
dftrain_filtrado = dftrainall_est[['GSPC.ReturnSuavizado','MSTF.ReturnSuavizado','BRK.B.ReturnSuavizado','VIX.ReturnSuavizado']]
dftest_filtrado = dftestall[['GSPC.ReturnSuavizado','MSTF.ReturnSuavizado','BRK.B.ReturnSuavizado','VIX.ReturnSuavizado']]
In [32]:
# Configuramos la función de pronóstico de forecaster con la mejor combinación obteida de hiperparámetros.
# ==============================================================================
### Aquí debe ir tu codigo, 
regressor = RandomForestRegressor(max_depth= 10, n_estimators=500 , random_state=123)
### Hasta aquí modificas

forecaster2 = ForecasterAutoreg(
                regressor = regressor,
                lags      = 5
             )

forecaster2.fit(y = dftrain_filtrado.iloc[:,0], exog = dftrain_filtrado.iloc[:,1:])
In [33]:
# Pronóstico, en este caso elegimos una ventana de 20 días.
# ==============================================================================
pred_filtrado = forecaster2.predict_interval(steps=steps, exog=dftest_filtrado.iloc[:,1:],interval = [1, 99],n_boot   = 500)
In [59]:
grafico_intervalos(dftest_filtrado.iloc[:,0],pred_filtrado)
desempeño(dftest_filtrado,pred_filtrado)

EL MSE del modelo en test es de: 0.14504318029190116
EL RMSE del modelo en test es de: 0.3808453495736835
EL R2 del modelo en test es de: 0.8906032684093637
EL MAE del modelo en test es de: 0.2971541961693256

Al evaluar las diferencias en las métricas encontramos que el R2 y el mse no cambian significativamente con el ajuste de las variables, por lo que continuaremos con la experimentación para reducir al máximo la cantidad de variables a utilizar.

En el siguiente experimento, quitaremos el VIX. Al ser el VIX un índice de volatilidad, en tiempo real, que se emplea como referente para cuantificar las expectativas del mercado en relación con la volatilidad implícita del S&P 500 (SPX) durante los siguientes 30 días, tener sus valores futuros implica realizar un proceso similar al que estamos haciendo para predecir el S&P500, por lo cual analizaremos el comportamiento de este modelo sin este índice para evaluar sus resultados:

Experimento # 3: Forecaster Autoregresor con variables exógenas filtradas (MSFT y BRK.B)¶

In [35]:
dftrain_sinvix = dftrain[['GSPC.ReturnSuavizado','MSTF.ReturnSuavizado','BRK.B.ReturnSuavizado']]
dftest_sinvix = dftest[['GSPC.ReturnSuavizado','MSTF.ReturnSuavizado','BRK.B.ReturnSuavizado']]
dftestall_sinvix = dftestall[['GSPC.ReturnSuavizado','MSTF.ReturnSuavizado','BRK.B.ReturnSuavizado']]
In [36]:
# Configuramos la función de pronóstico de forecaster con la mejor combinación obteida de hiperparámetros.
# ==============================================================================
### Aquí debe ir tu codigo, 
regressor = RandomForestRegressor(max_depth= 10, n_estimators=500 , random_state=123)
### Hasta aquí modificas

forecaster_sinvix = ForecasterAutoreg(
                regressor = regressor,
                lags      = 5
             )

forecaster_sinvix.fit(y = dftrain_sinvix.iloc[:,0], exog = dftrain_sinvix.iloc[:,1:])
In [37]:
# Pronóstico, en este caso elegimos una ventana de 20 días.
# ==============================================================================
pred_sinvix = forecaster_sinvix.predict_interval(steps=steps, exog=dftest_sinvix.iloc[:,1:],interval = [1, 99],n_boot   = 500)
In [60]:
grafico_intervalos(dftest_sinvix.iloc[:,0],pred_sinvix)
desempeño(dftest_sinvix,pred_sinvix)

EL MSE del modelo en test es de: 0.22132873517154628
EL RMSE del modelo en test es de: 0.47045588015407597
EL R2 del modelo en test es de: 0.8342890621228852
EL MAE del modelo en test es de: 0.3631235752483951

Aunque si perdemos desempeño, encontramos que la predicción del VIX resulta tan compleja como la predicción del S&P500, razón por la cual se contempla omitir esta variable.

Sin embargo, antes de tomar la decisión veamos qué tanto podría aportar al modelo únicamente el VIX:

Experimento # 4: Forecaster Autoregresor con variables exógenas filtradas (VIX)¶

In [39]:
dftrain_filtrado2 = dftrainall_est[['GSPC.ReturnSuavizado','VIX.ReturnSuavizado']]
dftest_filtrado2 = dftestall[['GSPC.ReturnSuavizado','VIX.ReturnSuavizado']]
In [40]:
# Configuramos la función de pronóstico de forecaster con la mejor combinación obteida de hiperparámetros.
# ==============================================================================
### Aquí debe ir tu codigo, 
regressor = RandomForestRegressor(max_depth= 10, n_estimators=500 , random_state=123)
### Hasta aquí modificas

forecaster_vix = ForecasterAutoreg(
                regressor = regressor,
                lags      = 5
             )

forecaster_vix.fit(y = dftrain_filtrado2.iloc[:,0], exog = dftrain_filtrado2.iloc[:,1:])
In [41]:
# Pronóstico, en este caso elegimos una ventana de 20 días.
# ==============================================================================
pred_filtrado2 = forecaster_vix.predict_interval(steps=steps, exog=dftest_filtrado2.iloc[:,1:],interval = [1, 99],n_boot   = 500)
In [61]:
grafico_intervalos(dftest_filtrado2.iloc[:,0],pred_filtrado2)
desempeño(dftest_filtrado2,pred_filtrado2)

EL MSE del modelo en test es de: 0.7542956626596483
EL RMSE del modelo en test es de: 0.868501964683816
EL R2 del modelo en test es de: 0.4310833506139946
EL MAE del modelo en test es de: 0.703442469348462

Este resultado nos muestra que si bien el VIX es importante, nos arroja una mayor precisión agregar los predictores de las acciones de MSTF y BRK.B, por lo cual la decisión final de los analistas será trabajar en dos modelos adicionales para predecir el comportamiento de estas dos acciones y omitir los datos del VIX como variable de entrada para el modelo

Resultados finales¶

Corramos el modelo con el 100% de los datos de entrenamiento ahora:

In [43]:
dftrainall_sinvix = pd.concat([dftrain_sinvix,dftest_sinvix])
dftestall_sinvix = dftestall[['GSPC.ReturnSuavizado','MSTF.ReturnSuavizado','BRK.B.ReturnSuavizado']]
In [44]:
# Configuramos la función de pronóstico de forecaster con la mejor combinación obteida de hiperparámetros.
# ==============================================================================
### Aquí debe ir tu codigo, 
regressor = RandomForestRegressor(max_depth= 10, n_estimators=500 , random_state=123)
### Hasta aquí modificas

forecaster_sinvix = ForecasterAutoreg(
                regressor = regressor,
                lags      = 5
             )

forecaster_sinvix.fit(y = dftrainall_sinvix.iloc[:,0], exog = dftrainall_sinvix.iloc[:,1:])
In [45]:
# Pronóstico, en este caso elegimos una ventana de 20 días.
# ==============================================================================
pred_sinvixfinal = forecaster_sinvix.predict_interval(steps=steps, exog=dftestall_sinvix.iloc[:,1:],interval = [1, 99],n_boot   = 500)
In [62]:
grafico_intervalos(dftestall_sinvix.iloc[:,0],pred_sinvixfinal)
desempeño(dftestall_sinvix,pred_sinvixfinal)

EL MSE del modelo en test es de: 0.18220597162151003
EL RMSE del modelo en test es de: 0.4268559143569526
EL R2 del modelo en test es de: 0.862573767814698
EL MAE del modelo en test es de: 0.37095388952853675

Tomemos entonces las predicciones obtenidas con el modelo sin el VIX para encontrar los resultados del modelo en retornos, recordemos que escalamos los datos y debemos volverlos a sus dimensiones originales:

In [47]:
# Construimos una matriz igual a la del conjunto de test para cada uno de los 3 intervalos de predicción
# ==============================================================================
#Matriz para predicciones
pred_filtrado_return = pd.DataFrame(0, index=pred_sinvixfinal.index, columns=dftestall.columns)
pred_filtrado_return['GSPC.ReturnSuavizado'] = pred_sinvixfinal['pred']
#Matriz para intervalo inferior
predlow_filtrado_return = pd.DataFrame(0, index=pred_sinvixfinal.index, columns=dftestall.columns)
predlow_filtrado_return['GSPC.ReturnSuavizado'] = pred_sinvixfinal['lower_bound']
#Matriz para intervalo superior
predup_filtrado_return = pd.DataFrame(0, index=pred_sinvixfinal.index, columns=dftestall.columns)
predup_filtrado_return['GSPC.ReturnSuavizado'] = pred_sinvixfinal['upper_bound']


#Transformamos los datos a la representación original.
pred_returns1 = pd.DataFrame(scaler.inverse_transform(pred_filtrado_return, copy=None))[0]
predlow_returns1 = pd.DataFrame(scaler.inverse_transform(predlow_filtrado_return, copy=None))[0]
predup_returns1 = pd.DataFrame(scaler.inverse_transform(predup_filtrado_return, copy=None))[0]

#Construimos un dataframe con las predicciones en la representación original
pred_final_1 = pd.concat([pred_returns1,predlow_returns1,predup_returns1],axis=1)
pred_final_1.columns = ["pred","lower_bound","upper_bound"]
pred_final_1.index = pred_sinvixfinal.index
In [48]:
# Importamos la información sin estandarizar de los retornos que tenemos en la predicción
# ==============================================================================
return_real = pd.DataFrame(df2.iloc[2792:2812,0])
return_real.head()
Out[48]:
GSPC.ReturnSuavizado
Date
2023-02-01 1.039806
2023-02-02 1.459238
2023-02-03 -1.040859
2023-02-06 -0.615939
2023-02-07 1.279036
In [63]:
# Calculamos las métricas
# ==============================================================================
grafico_intervalos(return_real.iloc[:,0],pred_final_1)
desempeño(return_real,pred_final_1)

EL MSE del modelo en test es de: 0.14076582767600526
EL RMSE del modelo en test es de: 0.3751877232479832
EL R2 del modelo en test es de: 0.8625737678146981
EL MAE del modelo en test es de: 0.32605228265810177
In [50]:
dftestall['GSPC.ReturnSuavizado'].max(),dftestall['GSPC.ReturnSuavizado'].min()
Out[50]:
(1.5971519617016467, -2.3663176740723677)

Finalmente, vamos a entrenar el modelo con todos los datos disponibles para futuras predicciones:

In [51]:
df_completo = pd.concat([dftrainall_sinvix,dftestall_sinvix])
In [52]:
# Configuramos la función de pronóstico de forecaster con la mejor combinación obteida de hiperparámetros.
# ==============================================================================
### Aquí debe ir tu codigo, 
regressor = RandomForestRegressor(max_depth= 10, n_estimators=500 , random_state=123)
### Hasta aquí modificas

forecaster_final = ForecasterAutoreg(
                regressor = regressor,
                lags      = 5
             )

forecaster_final.fit(y = df_completo.iloc[:,0], exog = df_completo.iloc[:,1:])

Exportamos el modelo final y el escalador¶

In [53]:
import joblib
joblib.dump(scaler, 'scaler.gz')
Out[53]:
['scaler.gz']
In [54]:
# Guardar modelo
#=====================================================================
save_forecaster(forecaster_final, file_name='forecaster.py', verbose=False)
In [55]:
load_forecaster('forecaster.py')

================= 
ForecasterAutoreg 
================= 
Regressor: RandomForestRegressor(max_depth=10, n_estimators=500, random_state=123) 
Lags: [1 2 3 4 5] 
Transformer for y: None 
Transformer for exog: None 
Window size: 5 
Weight function included: False 
Exogenous included: True 
Type of exogenous variable: <class 'pandas.core.frame.DataFrame'> 
Exogenous variables names: ['MSTF.ReturnSuavizado', 'BRK.B.ReturnSuavizado'] 
Training range: [Timestamp('2012-05-21 00:00:00'), Timestamp('2023-02-28 00:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: B 
Regressor parameters: {'bootstrap': True, 'ccp_alpha': 0.0, 'criterion': 'squared_error', 'max_depth': 10, 'max_features': 1.0, 'max_leaf_nodes': None, 'max_samples': None, 'min_impurity_decrease': 0.0, 'min_samples_leaf': 1, 'min_samples_split': 2, 'min_weight_fraction_leaf': 0.0, 'n_estimators': 500, 'n_jobs': None, 'oob_score': False, 'random_state': 123, 'verbose': 0, 'warm_start': False} 
Creation date: 2023-05-07 16:26:00 
Last fit date: 2023-05-07 16:26:09 
Skforecast version: 0.7.0 
Python version: 3.9.13 
Forecaster id: None
Out[55]:
================= 
ForecasterAutoreg 
================= 
Regressor: RandomForestRegressor(max_depth=10, n_estimators=500, random_state=123) 
Lags: [1 2 3 4 5] 
Transformer for y: None 
Transformer for exog: None 
Window size: 5 
Weight function included: False 
Exogenous included: True 
Type of exogenous variable: <class 'pandas.core.frame.DataFrame'> 
Exogenous variables names: ['MSTF.ReturnSuavizado', 'BRK.B.ReturnSuavizado'] 
Training range: [Timestamp('2012-05-21 00:00:00'), Timestamp('2023-02-28 00:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: B 
Regressor parameters: {'bootstrap': True, 'ccp_alpha': 0.0, 'criterion': 'squared_error', 'max_depth': 10, 'max_features': 1.0, 'max_leaf_nodes': None, 'max_samples': None, 'min_impurity_decrease': 0.0, 'min_samples_leaf': 1, 'min_samples_split': 2, 'min_weight_fraction_leaf': 0.0, 'n_estimators': 500, 'n_jobs': None, 'oob_score': False, 'random_state': 123, 'verbose': 0, 'warm_start': False} 
Creation date: 2023-05-07 16:26:00 
Last fit date: 2023-05-07 16:26:09 
Skforecast version: 0.7.0 
Python version: 3.9.13 
Forecaster id: None